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We investigate the energy transfer dynamics in a donor-acceptor model by developing a time-local master 
equation technique based on a variational transformation of the underlying Hamiltonian. The variational 
transformation allows a minimisation of the Hamiltonian perturbation term dependent on the system pa- 
rameters, and consequently results in a versatile master equation valid over a range of system-bath coupling 
strengths, temperatures, and environmental spectral densities. While our formalism reduces to the well-known 
Redfield, Forster and polaron forms in the appropriate limits, in general it is not equivalent to perturbing in 
either the system-environment or donor-acceptor coupling strengths, and hence can provide reliable results 
between these limits as well. Moreover, we show how to include the effects of both environmental correlations 
and non-equilibrium preparations within the formalism. 



I. INTRODUCTION 

The transfer of an electronic excitation is a ubiqui- 
tous process in physics, chemistry, and biology. Typi- 
cally, coupling to the radiation field creates an excitation 
in one location (the donor), and through the exchange 
of a virtual photon, the excitation is passed to another 
(the acceptor).^ One of the many challenges in theoreti- 
cally modelling such a process lies in correctly accounting 
for the influence of the environment surrounding the two 
(or more) sites, which determines whether the transfer 
is of a coherent or an incoherent nature. A common ap- 
proach has been to assume that the energy transfer cou- 
pling strength is weak, and thus perform perturbation 
theory in this parameter. The resulting Forster-Dexter 
rate equations then describe purely incoherent energy 
transfer, ^'"^ and are applicable in a wide variety of situa- 
tionsj^i^ However, recent experimental results providing 
evidence for quantum coherent transfer in a number of 
systems^"— have necessitated the development of theo- 
retical techniques beyond this approximation. While ex- 
tending Forster theory to account for coherence within 
multi-site donors or acceptors is possible,J^^"^^ describing 
the dynamical evolution of coherences between the donor 
and acceptor generally requires an alternative starting 
point. 

One such approach is to treat the system-environment 
coupling perturbatively, while accounting for the elec- 
tronic coupling between the sites to all orders. Work- 
ing with the resulting Redfield or Lindblad equations 
is vastly simpler than simulating the full dynamics)^ 
and has thus recently been put to great effect to anal- 
yse the interplay of dissipation and quantum coher- 



ence in the energy transfer dynamics of complex, many- 
site systems j^"""2P, However, by their very nature, treat- 
ments based upon a weak system-environment coupling 
assumption are often inadequate to describe systems 
strongly coupled to their environments and/or at high 



temperatures 
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In this situation, there exists a num- 



^) Electronic mail: dara.mccutcheon@ucl.a c.ukl 
'"'Electronic mail: ,a.nazir@imperial. ac. uk| 



ber of (non-perturbative) numerical approaches which al- 
low one to explore the full range of system-bath cou- 
pling and electronic transfer strengths, subject to cer- 
tain method-specific constraints. Examples include den- 
sity matrix renormalisation group;^ numerical renormal- 
isation group, "^^ and path integral^ii^S methods. Addi- 
tionally, for specific forms (but not strengths) of system- 
environment coupling, numerically exact results can be 
obtained through the hierarchical equations of motion 
technique)^!^"— which has recently been shown to be 
consistent with the path integral formalism4^ 

These methods have clear advantages in that their 
range of applicability is generally less restricted than 
those of the aforementioned perturbative approaches. 
However, given the physical insight perturbative tech- 
niques can impart, and the relative simplicity and effi- 
ciency with which they can be implemented, it remains 
of particular importance to develop approximation meth- 
ods which allow one to probe as wide a range of pa- 
rameter regimes as possible. For example, aspects of 
energy transfer dynamics have recently been explored 
using a time-local polaron-transform master equation 
technique^"— which corresponds to perturbation theory 
in an environment-dressed electronic coupling ternij^"— 
For particular forms of the system-environment coupling, 
the polaron master equation reduces to Redfield the- 
ory in the weak-coupling limit, but also remains valid 
for far greater temperatures and system-bath coupling 
strengths, ^Si^i^ and can thus also capture the Forster 
limit. This allows, in particular, for a consistent ex- 
ploration of the crossover from coherent to incoherent 
dynamics in energy transfer processes4^^ The polaron 



method, however, suffers itself from a restriction to rela- 
tively small electronic coupling strengths (compared to 
typical frequencies in the environment), and, as men- 
tioned, can only interpolate between the Redfield and 
Forster limits for certain forms of system-environment 
coupling. 

In this work we go beyond the polaron approach and 
present a master equation technique which combines a 
variationally-optimised unitary transformation^^— with 
the time-convolutionless projection operator formal- 
ism J^ Modelling the environment as a bath of harmonic 
oscillators, the variational transformation displaces each 
mode according to the position of the excitation, and 
by an amount defined to minimise the resulting free en- 
ergy of the interaction terms between the donor-acceptor 
pair and the environment. We contrast the resulting 
time-local master equation with the Redfield and polaron 
forms, and show that the variational approach is supe- 
rior to both. Specifically, we show that while the varia- 
tional master equation reproduces both the Redfield and 
polaron equations in the appropriate limits, it can also 
give qualitatively reliable results outside these parameter 
regimes. The variational master equation can therefore 
be used to explore energy transfer dynamics for a range of 
system-environment coupling forms, strengths, and envi- 
ronmental temperatures. 

The paper is organised as follows: In Section [11] we 
introduce the donor-acceptor model and the variational 
transformation. We then outline the derivation of the 
time-local master equation in Section IIIII In Section IIVI 
we use this master equation to explore the influence 
of both super-Ohmic and Ohmic environments on the 
donor- acceptor energy transfer dynamics. We also con- 
sider here the role of bath relaxation effects, and how 
they influence the rate of energy transfer. A brief sum- 
mary is presented in Section |Vl while Appendix |X] gives 
further details of the master equation derivation, and Ap- 
pendix |B] extends the formalism to consider correlated 
environmental fluctuationsi^i^'^i^'^".^ 



II. MODEL AND VARIATIONAL TRANSFORMATION 



We consider a donor-acceptor pair {j — 1, 2) each site 
of which we model as a two-level system with ground 
state \G),- and excited state \X)-, split by an energy hej. 
The pair interact via an electronic coupling of strength V 
which transfers an excitation from one site to the other. 
To model the dephasing and dissipative effects of the en- 
vironment, we linearly couple each excited state to a bath 
of harmonic oscillators i^^^iS Although not the primary 
focus of this work, the formalism to be presented below 
can also been extended to include the effects of correlated 
fluctuations between the sites.'**'** This extension is out- 
lined in Appendix [BJ In the absence of correlations, the 



total Hamiltonian reads (we set h = 1) 

i/ = ^ 6, \X)^{X\ + V{ \XG){GX\ + \GX){XG\ ) 

+ ^ \X)^{X\J29kAbh+h^j) + HB, (1) 

i=l,2 k 

where \X)^{X\ = \X){X\ ® {\X){X\ + \G){G\) and 
\X)^{X\ = i\X){X\ + \G){G\) (g> \X){X\ determine the 
excited state populations of the two sites, gkj is the 
coupling constant of site j to mode k of bath j, de- 
scribed by creation (annihilation) operators bj, (bk.j)- 
The bath Hamiltonian is Hb — J2jJ2k'^k.jbk jbkj, 
with frequencies Ukj- Eq. ([T]) generates excitation dy- 
namics within the three decoupled subspaces, spanned 
by {\GG) ,{\XG) ,\GX)},\XX)}. We are interested 
here in the single-excitation subspace (spanned by 
{\XG),\GX)}) and therefore set \XG) = |1) and 
\GX) = 1 2), and write the Hamiltonian governing dy- 
namics within this subspace as 



^STIR — -CTj: 



V(J^ + He 



+ E \jMj2!^''M,,+bk,,). (2) 

j=1.2 k 

Here, e = ei — 62 is the energetic difference between the 
two sites, the Pauli operators are defined in the basis 
az = |1)(1| — |2)(2|, and a term proportional to the iden- 
tity has been neglected. 

A standard way to proceed from this point might be 
to treat the system-bath interaction term in Eq. ([2]) as a 
perturbation, resulting in a master equation of Lindblad 
or Redfield typeJ^ Alternatively, a polaron transforma- 
tion could first be applied to Eq. ([2]), allowing one to 
identify a new interaction term, and result ing in a mas- 
ter equation of the form explored in Refs. I44l448l . Here, 
we instead employ a method originally developed by Sil- 
bey and Harris and apply a variational transformation 
to i?suB-— "^ As in the polaron case, the transformation 
displaces those modes in the bath coupled to the site 
possessing the excitation. However, unlike the full po- 
laron transformation, we allow freedom within the vari- 
ational transformation to attempt to optimise these dis- 
placements for each mode. This will be achieved through 
free-energy minimisation arguments 1^""— 

The transformed Hamiltonian is defined by Ht = 
e'^^suBe"'^, with 






'fcj 



bk,3 



(3) 



This results in Ht = Hq + Hj, with free Hamiltonian 

Ho = i(e + i?i-i?2)a, + yflfT,+i/s + i(i?i+i?2)l, (4) 

where Rj = J2k hd^kjih,] - '^9k,j), and the renor- 
malised electronic coupling strength Vr will be defined 



below. The interaction Hamiltonian now contains two 
terms, Hi = -ffLiNEAR + -f^DiSPLACED, with the first given 
given by 

-ffuNEAR = ^ |j)01 ^{gk,j - /fc,j)(^Ij + bk,j), (5) 
3 k 

which has the same form as the perturbative term in 
Redfield theory, but with the coupling to each mode now 
reduced in strength. Notice that if we were to set fk,j = 
gu.j for all modes, as in the full polaron transformation, 
-f^LiNEAR would vanish as expected, though in general 
this is not the case. The second term in Hj has the form 
of the perturbation used in deriving a polaron master 
equation, and is given by 



-f^DISPLACED — V{(Tj;Bx 



ay By) 



(6) 



written in terms of the Hermitian combinations B^ — 
(1/2) (B+ + B_ - 2B) and By = {i/2){B+ - B_) of the 
bath displacement operators operators B± — -B±,ii?^.2, 
with 



B±,,=exp\±J2l^{bl^~b,/, 



(7) 



Again, it is important to note that only in the particular 
case of fk,j = gk,j (for all modes) does Eq. ^ become 
identical to the polaron form. 

Going back to the free Hamiltonian, we see that the 
term generating the energy transfer in Eq. Q now ap- 
pears with a renormalised strength; Vr = VB, with 
B = Ti{B±pb) being the expectation value of B± with 
respect to the bath state. Taking a thermal equilibrium 
state pb ~ e~'^^^/(Tre~^^^) at inverse temperature 
/3 ~ l/ksT gives explicitly 



B 



■ exp 



EE 



f? 



kj 



coth 



/ P^fcj \ 



(8) 



If we now assume that the baths coupled to the donor 
and acceptor are identical, this allows us to set gk^i = 
gk,2 = gk, 0Jk,i = Wfc,2 = Wfc, /fc,i = fk^2 — fk, and i?i = 
R2 = R = ^k fk'^k^ifk ~ '^gk)- The free Hamiltonian 
now takes on the simpler form 



Hn = -a. 



VRa,.,+Rt + HB, 



while the renormalisation factor becomes 



B = exp 



[-E 



k ^fc 



I coth(/3Ljfc/2) 



(9) 



(10) 



A. Free energy minimisation 

Currently, the parameters {/fc} appearing in our trans- 
formed Hamiltonian are undetermined. We note that set- 
ting fk — corresponds to having performed no transfor- 
mation on the Hamiltonian, and as such -ffoiSPLACED = 



while -ffLiNEAR remains finite, and the Hamiltonian takes 
on its original spin-boson form.^^ On the other hand, as 
remarked earlier, by setting fk ~ gk for all k one finds 
that i^LiNEAR = while i^oisPLACED remains finite. In 
this case, the Hamiltonian corresponds to that of a po- 
laron transformed spin-boson model4i^— i^i^ Our aim 
in this work is to attempt to derive a second-order (time- 
local) master equation valid over as large a range of pa- 
rameter space as possible. We achieve this by trying to 
find an interaction Hamiltonian which remains small over 
the greatest range of parameters. Therefore, rather than 
setting /it = or /fe = (7fc, we instead determine the set 
{fk} by free energy minimisation arguments. 

To proceed, we compute the Feynman-Bogoliuobov up- 
per bound on the free energyj^i^ given by 



Af 



ln(Tr{e-^^n) + {Hi)h,. + 0{{HJ)h^X (H) 



and related to the true free energy, A, via the inequality 
Ab > A^ By construction, the second term appearing 
in Eq. (jlip is equal to zero. The variational transforma- 
tion does not affect the value of the free energy of the 
complete system; in order to minimise the contribution 
from Hi, we therefore neglect the third term in Eq. (fTTj) 
and minimise (or maximise in magnitude) the first term 
with respect to the variational parameters. We find 



^B = i?-^ln[2cosh(/3r//2)], 



(12) 



with rj = -^/e^ + 4V^, and minimisation with respect to 
fk gives fk = gkF{ujk), with 



F{ujk) 



1 + ^ tanh (/J77/2) coth (/3a;fe/2) 
r]ujk 



(13) 



Thus, in general, the interaction Hamiltonian contains 
contributions from both -//linear and ii/oisPLACED • As 
we shall see below, in deriving a master equation utilising 
Hj as a perturbation, we therefore have terms arising 
from both i^LiNEAR and i/oisPLACED, as well as from 
their product. 

Introducing the spectral density (assumed to be the 
same at each site) as J{uj) = J2k \9k\'^S{^ ~ ^k) allows 
Eq. (fTU|) to be written in integral form as 



B 



exp 



duj 



J(c.) 



F{ujfcoth{(3uj/2) . (14) 



Since Vr — VB, and B is itself a function of Vr, the 
renormalised coupling strength must be solved for self- 
consistently. 

For a sufficiently large bath of oscillators the spec- 
tral density is typically taken to be a smooth func- 
tion of w, with polynomial-like behaviour in the small 
uj limit; J(w) ^ cj^ as a; — > 0. Those spectral densi- 
ties with s < 1, s = 1, and s > 1 are referred to as 
sub-Ohmic, Ohmic, and super-Ohmic, respectively)^ Of 
particular significance is the value of the renormalised 



interaction strength, Vr, found for an Ohmic environ- 
ment. For the full polaron transformation, the integral 
in Eq. (|14p suffers from a well-known infra-red diver- 
gence,^^'^^ which can be seen here by setting F{uj) — 1 
and taking J{uj) ~ w. In this case, Vr — > in the po- 
laron transformation regardless of the temperature or the 
strength of the system-bath coupling. For an Ohmic en- 
vironment, we are therefore forced to conclude that only 
incoherent dynamics can be captured by the full polaron 
transformation when it is used in conjunction with a 
time-local master equation approach. Importantly, this 
is not the case for the variational theory as the function 
F{u}) need not be equal to 1, and hence both coherent 
and incoherent dynamics can be captured. 



III. MASTER EQUATION 

Having performed the variational transformation on 
our Hamiltonian to identify the perturbation term Hj, 
we now employ the time-convolutionless projection op- 
erator methodiS. to derive a master equation governing 
the donor-acceptor energy transfer dynamics under the 
influence of the surrounding environment. This tech- 
nique utilises a projection operator, 7-", defined to sat- 
isfy Px = TrB(x) ^ p-R = p <S) PR, where x is the den- 
sity operator of the combined system-plus-environment 
state, p is the reduced density operator of the system 
(the donor-acceptor pair in our case), and pR is a ref- 
erence state of the environment, which can in principle 
be chosen arbitrarily. Using the projection operator, an 
exact time-local master equation for Vx can be derived, 
which involves Vx itself and, in general, the complemen- 
tary projection, Qx(0) = (-^~'^)x(0)i of the initial state 
x(0). Importantly, in the present case, since our Hamil- 
tonian has been transformed into the variational frame, 
we obtain a master equation involving the variationally- 
transformed density operator, xt = e'^x^~'~^- Truncat- 
ing the exact expression to second order in Hj, our in- 
teraction picture master equation reads 

^ = Ti-B[mVxT{t)]+TrB[mQxTm, (15) 
where the homogeneous contribution is given by 



TrB[ICit)PxT{t)] = - / dsTTB[Hi{t),[Hi{s),pT{t)pR]], 
Jo 

(16) 

while the initial-state-dependent inhomogeneous term is 
TvBimQxTm = -tTr b[H I {t),QxTm 



dsTi-B[Hj{t),[Hi{s),QxTm], 

(17) 



with tildes indicating operators in the interaction pic- 
ture, Hj{t) = eyip[iHot]Hiexp[—iHot], and prit) — 
Tr B{xT{t)) being the reduced density operator in the 



variationally-transformed frame. Note that in deriving 
Eq. ([T5|) we have utihsed the fact that TTB{Hi{t)pR) = 0. 
Derivation of the final form of the master equation 
now proceeds in the usual way, and is somewhat lengthy 
given the various terms present in Hj . We therefore leave 
the details for Appendix \X\ and mention only the salient 
points here. From Eqs. ([TS]) and ([TT]) it is evident that the 
master equation will contain two-time correlation func- 
tions arising from the product of Hj{t) = i?LiNEAR(^) + 
-ffDiSPLACED(0- We write Hj = J2i=i ^i ® Bi, with 
system operators Ai — |1)(1|, A2 — |2)(2|, A3 — a^, and 
yl4 = cTj^, and bath operators Bi = I]fc(3fc~/fc)(^L+&fc,i) 
for i = 1,2, while B3 = VB^ and B4 = VBy. The cor- 
relation functions appearing in the master equation can 
then be written as 



A„(r)=TrB(B,(T)Bj(0)pR), 



(18) 



where we choose a thermal equilibrium bath reference 
state, PR = e"'^^«/(Tre~'^^«). 

For i = 1,2 these correlation functions have a similar 
form to those found in Redfield theory: 



Aii(t) = A22(r) 



da;J(w)(l -F{u)f 



X (cos(a;r) coth(/3aj/2) — i sin(ajT)) , 

(19) 

while Ai2(r) = A2i(t) — 0. Note that in the polaron 
limit {F{uj) -^ 1) these correlation functions vanish. For 
i — 3,4, we obtain correlation functions similar in form 
to those of a polaron master equation)^"— 

A33(t) - (V-^/2)(e*M + e-^(-) - 2), (20) 



A44(r) = (F|/2)(c^(^ 



-<t>ir 



(21) 



with 



f° 
(r) =2 / 

Jo 



du! 



J(u;) 



F{^)' 



X (cos(wt) coth{f3uj/2) — i sin(a;r)) 



(22) 



while A34(t) = A43(t) = 0. In the weak-coupling (Red- 
field) limit F{uj) — i> 0, and A33(t) and A44(t) then dis- 
appear. The last type of correlation function is unique 
to the variational theory and arises from the product of 
^LINEAR and Hdisplaced- We have Ai4(t) = A42(t) = 
iAcir), and A24(t) = A4i(r) = —iAcir), where 

Ac(r) = Vn H du:^^F{Lo)il - F(c.)) 



X (sin(cjT) coth(;9w/2) + i cos(wt)) , (23) 

and A3i(t) = Ai3(t) = A32(t) = A23(t) = 0. These 
correlation functions are most important in the interme- 
diate regime, where neither the full polaron displacement 
{F{uj) — 1) nor zero displacement {F{uj) — 0) are appro- 
priate. 



The inhomogeneous term in Eq. p5|) depends on the 
quantity Qxt(O)- Assuming a separable initial state cor- 
responding to an excitation on the donor, p{0) = |1)(1|, 
with the environment in the state pb{0), wc find 



Qxt(O) = (1 - P)(e« |1)(1| ® PB(0)e-^) 
= |1)(1|®(ptb(0)-pr), 



(24) 



where ptb(O) = -B+.iPb(0)-B_.i is the variationally- 
transformed bath initial state. Thus, we see that the 
inhomogeneous term vanishes if we assume our initial 
environmental state to be displaced such that /9b(0) = 
B-^iPY\B+,i, since then ptb(O) = pr. We might, how- 
ever, wish to consider a case where the excitation enters 
the system on a short enough timescale such that the 
bath does not have time to displace in response to it. The 
correct initial state would then be pb{0) = PR, such that 
Ptb(O) = B+.iPrB-,!. We then find Qxt(O) 7^ 0, and so 
the inhomogeneous term remains. Details of the first and 
second-order inhomogeneous contributions to the master 
equation in this case are presented in Appendix [Xl 
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FIG. 1. Excitation dynamics for various reorganisation ener- 
gies and cut-off frequencies, and for coupling to a super-Ohmic 
environment. Calculations using the variational master equa- 
tion (black solid curves), Redfield theory (red dashed curves), 
and polaron theory (blue dotted curves) are shown. Parame- 
ters: e = too cm"\ V = too cm"^ and T = 300 K. 



IV. ENERGY TRANSFER DYNAMICS 

Let us now use the variational master equation we 
have derived to investigate the dynamics of an excitation 
in the donor-acceptor system. To do so, we calculate 
the population in state |1) (the donor) as a function of 
time, pii(t) = (1/2)(1 + a,{t)), where a, = Ti{a,x{t))- 
Generally, care must be taken when calculating expec- 
tation values from a master equation derived in a trans- 
formed frame. However, since e aze~'^ = (Tz, we find 
a,{t) = Tv{a,xit))= Tr(a,e-GxT(i)e«) = Tr{a,XT{t)). 
Thus, the population dynamics in which we are inter- 
ested is unaffected by transformation into the variational 
representation, and we may use our master equation for 
pxit) directly to compute pii(t). 

In the following sections we shall compare the dynam- 
ics generated by the present variational theory to that 
found using both the Redfield and polaron approaches. 
Redfield theory corresponds to second-order perturba- 
tion in the original system-bath interaction term, before 
any transformation has been applied to the Hamiltonian. 
Within the present formalism it can therefore be recov- 
ered by setting F{uj) = 0. Typically, it is further assumed 
in this limit that the bath correlation functions decay on 
a timescale much faster than that of the excitation dy- 
namics, which allows the integration limits in Eq. (J16p 
(and also in Eq. (jA3[) ) to be taken to infinity. Dynamics 
referred to as Redfield in the subsequent sections is cal- 
culated in this way. On the other hand, a polaron master 
equation is obtained through performing the full trans- 
formation on the original Hamiltonian. In our formalism, 
this corresponds to setting F{uj) = 1. Making no further 
assumptions, our theory then reduces to that presented 
in Refs. |45| and |46|. Dynamics referred to as polaron in 
the following is calculated in this manner. 



A. Super-Ohmic environments 

To begin our analysis we shall first consider a donor- 
acceptor pair that is coupled to a super-Ohmic environ- 
ment, as has been studied previously using the polaron 
formalism in Refs. |44| - |49| . This immediately allows us to 
investigate how the present variational master equation 
theory compares to both the Redfield and polaron forms, 
without having to worry about complications such as the 
polaron theory infra-red divergence, which would be an 
issue in the Ohmic case (see Section FlVBp . We take a 
spectral density of the form 



Jsi^) 



a3w3e""/"<= 



(25) 



where a^ captures the strength of the system-bath cou- 
pling and uJc is a phenomenological cut-off frequency. 
A spectral density of this form is typical in the solid- 
state, for example when describing coupling to acoustic 
phonons, see e.g. Refs. l69l - [7ll We also introduce the 
reorganisation energy, defined as 



A. = 



Js{io)uj dcj. 



(26) 



which constitutes a measure of the system-bath interac- 
tion that also accounts for the range of frequencies over 
which the bath can influence the system. For the super- 
Ohmic spectral density given above we find A3 = 2a3Uj^. 
In Fig. [T] we plot the excitation population dynamics 
using the variational theory (solid black curves). Red- 
field theory (red dashed curves), and the polaron the- 
ory (blue dotted curves) , for various values of the reor- 
ganisation energy, A3, and cut-off frequency, ujc- In the 
variational and polaron cases we have assumed that the 
environment is initially in a displaced thermal equilib- 
rium state, resulting in the absence of inhomogeneous 



terms, as discussed in Section IIIIl The effect of assum- 
ing a non-displaced initial bath state is discussed in Sec- 
tion IIVCI In plot (a) the reorganisation energy has a 
relatively small value of A3 = 20 cm~^, and the cut-off 
frequency is such that V juJc < 1- For these parameters 
we expect both the Redfield and polaron theories to be 
valid,^':54 and hence they yield almost identical results. 
We also see that the variational theory predicts essen- 
tially the same dynamics in this relatively undemanding 
weak-coupling regime. 

Plot (b) corresponds to the case where the reorgani- 
sation energy has been increased by a factor of 10. We 
would not expect Redfield theory to be justified in this 
regime, owing to the large value of A3, and indeed we see 
that it predicts a strong damping of the system dynamics, 
at odds with the other approaches. In fact, in this regime 
the variational minimisation condition corresponds ap- 
proximately to performing the full polaron displacement, 
i.e. F{ujk) ~ 1, and so the polaron and variational theo- 
ries give very similar results. This is to be expected (for 
a super-Ohmic bath), as for large A3 and small V/lOc the 
full polaron transformation provides a much smaller per- 
turbative term than the original system-bath coupling;^ 
and is therefore a superior representation in which to ap- 
proximate the system dynamics. 

For the parameters of plot (c), on the other hand, we 
now have V/ujc > 1, and the full polaron displacement is 
no longer appropriate. We therefore see that the polaron 
theory incorrectly predicts a strong damping of the sys- 
tem dynamics, even though the reorganisation energy is 
small. In fact, the variational minimisation condition cor- 
responds here to performing only a weak transformation 
on the Hamiltonian, i.e. F{ujk) ~ 0, and the variational 
and Redfield theories thus predict similar results. 

Plot (d) corresponds to none of the limiting cases dis- 
cussed above, and consequently all three theories pre- 
dict different dynamics. Here, we expect the variational 
approach to most accurately describe the true dynam- 
ics, since the minimisation condition allows for a mode- 
specific optimisation of the displacement of the bath, 
which the other theories do not. We note also that a 
recent comparison of the three theories considered here 
with a numerically-exact path integral approach con- 
firmed the superiority of the variational method over 
Redfield and polaron for super-Ohmic environments, al- 
beit for a different model systemi^^ 



B. Ohmic environments 

We now consider the excitation dynamics when the 
donor-acceptor pair are coupled to an Ohmic bath. 
Specifically, we take the overdamped Brownian oscilla- 
tor form 




FIG. 2. Excitation dynamics calculated using the variational 
theory for coupling to an Ohmic environment, with reorgan- 
isation energies as indicated (in units of cm~^). Parameters: 
e = 100 cm"\ V = 100 cm"\ ujc = 53 cm"\ and T = 300 K. 
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FIG. 3. Redfield and polaron dynamics for coupling to an 
Ohmic environment. Each curve corresponds to a different 
reorganisation energy, as indicated in Fig. (2] Parameters: 
e = 100 cm"\ V = 100 cm"\ Wc = 53 cm"\ and T = 300 K. 



Ji(c^) 



2 aiLOcLo 



(27) 



which behaves linearly in the limit cj — )• 0. From Eq. ([26]) 
we find Ai = ai , which is equivalent to the reorganisation 
energy used in Ref. |3J. 

In Fig. [5] we plot the excitation dynamics calculated 
from the variational theory for several values of the cou- 
pling Ai. As in the super-Ohmic case, we assume a dis- 
placed initial bath state resulting in the absence of in- 
homogeneous terms. As we might expect, for small Ai 
the dynamics shows pronounced oscillations with a de- 
caying envelope, and as Ai is increased, the oscillations 
become more strongly damped. Once Ai = 100 cm~^, 
however, the dynamics becomes entirely incoherent, and 
the steady state is reached on a timescale ~ 0.1 ps. In- 
terestingly, as the coupling strength is increased further, 
the excitation begins to take longer to relax towards the 
steady state, resulting in a decrease in the rate of ex- 
citation transfer between the two sites. This behaviour 
is consistent with results obtained using the hierarchal 
equations of motion technique.— 1^ 

By way of comparison, in Fig. |3] we plot excitation dy- 
namics for parameters identical to those in Fig.[5J though 
now calculated using Redfield and polaron theories. As 
in the super-Ohmic case, we see that for small reorgani- 
sation energies, the Redfield and variational theories pre- 
dict very similar dynamics. However, as the system- 
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FIG. 4. Dynamics of the populations of BChl 1 (black) and 
BChl 2 (orange) for parameters corresponding to the FMO 
complex, at r = 77 K (left) and T = 300 K (right). Pa- 
rameters: Ai = 35 cm~^, LOc = 106 cm^^, e = —120 cm~^, 
V ^ -87.7 cm-\ 



bath coupling strength becomes large, Redfield theory 
cannot capture the expected reduction of the transfer 
rate^^i^ and there appears to be a particular coupling 
strength after which the steady state is always reached 
on a timescale '^ 0.1 ps. 

For the Ohmic spectral density studied here, we see 
that polaron theory behaves quite oddly. As previously 
mentioned, a spectral density of this form presents a 
problem for a time-local master equation obtained in the 
polaron frame, since a complete renormalisation of the 
electronic transfer strength is predicted, Vr — ;■ 0. In this 
situation, the polaron formalism can capture only inco- 
herent dynamics. In fact, when Vr ~ (in either the 
polaron or variational theory), from Eq. (|A1[) we find 
a simple expression governing the time evolution of the 
population on site 1: 



dpiijt) 
dt 



= -At(e,t)pn(i) + «:(-£, i)(l - pii(i)), (28) 



where the rates determining transfer between the sites 
are given by2^ 



K(±e,i) = 2y2Re 



dre 



zLier 



('52g0(r)\ 1 (29) 



Thus, in the limit Vr ^ we see that the polaron theory 
reduces to Forster theory, though with time-dependent 
transfer rateS)^ which we expect to work well only in 
the strong system-environment coupling limit i^i^ This is 
confirmed by the unphysical behaviour we see in the po- 
laron plot for Ai = 1 cm~^. At larger coupling strengths, 
however, polaron theory does correctly predict a reduc- 
tion of the transfer rate. On comparison of Figs. [2]and|31 
it can be seen that the minimisation condition [Eq. ([T3| ] 
present in the variational theory picks out an optimised 
Hamiltonian transformation dependent upon the system- 
bath parameters, resulting here in qualitatively correct 
dynamical behaviour across the full range of coupling 
strengths. 

To test the variational theory in a definite physical 
context, in Fig. 2] we plot the energy transfer dynam- 
ics in a very simple model that consists of two strongly- 
coupled bacteriochlorophyll sites (BChl 1 and BChl 2) of 



the Fenna-Matthews-Olson (FMO) complexi^i^ We use 
our two-site Hamiltonian, Eq. ([1]), with system parame- 
ters taken from Ref. [7a (corresponding to e = —120 cm^^ 
and V = —87.7 cm~^, see also Ref. 77), and consider 
an Ohmic spectral density [Eq. (|27l) ] with reorganisation 
energy Ai = 35 cm^^ to match Ref. |40. By compari- 
son of our variational results at cryogenic temperature 
(T — 77 K) and physiological temperature (T = 300 K) 
with Figs. 2 and 3, respectively, in Ref. |40l it can be 
seen that the variational theory is in excellent agree- 
ment with the numerically-exact hierarchical approach 
on a timescale of the order ~ 100 fs, after which time we 
would expect discrepancies as population begins to leak 
into the other BChl sites not accounted for here 4^121"— 
Hence, while this may be a very simplified example, it 
serves to highlight the versatility of the variational ap- 
proach, and paves the way for a more detailed study of 
the FMO system. 



C. Bath relaxation effects 

In the previous sections we assumed a displaced ini- 
tial bath state when calculating dynamics using the vari- 
ational and polaron theories. This simplified the rel- 
evant master equations, since it meant that they con- 
tained no inhomogeneous terms. This simplification re- 
lies on a separation of timescales in the combined system- 
environment dynamics, i.e. the bath relaxation must 
be fast compared to the energy transfer dynamics be- 
tween the two sites. Thus, though the cut-off frequency 
dependence in the variational theory was seen to play 
an important role in the reduced dressing of the energy 
transfer interaction strength, the bath was still assumed 
to instantaneously adjust to its variationally displaced 
state. We now relax this assumption by the introduc- 
tion of the inhomogeneous terms, the presence of which 
results from a difference between the initial bath state 
and that taken to be the reference state in the projection 
operator, and therefore (approximately) account for the 
influence of environmental relaxation on the state of the 
system. Hence, we shall now consider a non-displaced 
initial bath state, pb(0) = pR, and investigate what dif- 
ference the inhomogeneous terms in our variational mas- 
ter equation may make, as has been studied previously 
in the polaron case4^i^'^ 

In Fig. [5] we plot the excitation dynamics calculated 
from the variational master equation for a super-Ohmic 
spectral density [Eq. ([25]) ] including (dashed curves) and 
excluding (solid curves) all first- and second-order inho- 
mogeneous terms. The main plot is for zero temperature, 
while the inset corresponds to T = 300 K. Though we 
find here that the effect of the inhomogeneous terms for 
a localised initial state is weak^ interestingly, at zero 
temperature we see that the relaxation of the bath ac- 
tually causes oscillations in the population dynamics to 
be suppressed. At T = 300 K, the situation is somewhat 
different, and the dynamics is predominately incoherent 
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FIG. 5. Excitation dynamics with (dashed) and without 
(solid) bath relaxation terms, for a super-Ohmic spectral 
density. The main plot corresponds to zero temperature, 
while the inset shows dynamics for T = 300 K. Parameters: 
A3 = 50 cm-\ LOc = 53 cm-\ e = 100 crcT^ , V = 100 cm'^ 



without inhomogeneous terms, whereas their inclusion 
seems to give rise to small amplitude oscillations.'*^ 

As a brief aside, at this point it is worth comparing 
the steady state reached at zero temperature in the vari- 
ational theory to the incorrect value predicted by the 

non-interacting blip approximation (NIBA) for a biased 
„Nir- 



systemi^iS pf/B^ = (l/2)(l-tanh(/3e/2)). At zero tem- 



perature, we find p^l^^ — 0, and so the NIBA predicts 
that all population is transferred to site 2, regardless of 
the reorganisation energy, or the relative values of V and 
e. From Fig. [S] we see that this is not the case in the 
variational theory, and we therefore expect it to better 
capture the corresponding dynamics in this regime. 



Energy Transfer Rate 

The effect of the bath relaxation terms can also be 
seen in the inter-site energy transfer rates, which we 
explore here for an Ohmic bath [Eq. (P7| ]. To calcu- 
late these rates, we fit the time domain dynamics to 
the solution of the simple classical rate equation pn = 
—K+pii + K_(l — pii), and determine the quantity k+, 
which characterises the rate at which excitation is passed 
from site 1 to site 2. Fitting quantum dynamics to such a 
rate equation is clearly dubious if there are many coher- 
ent oscillations present in the system. We shall therefore 
consider a small donor-acceptor interaction strength of 
y = 20 cm^^, for which we expect the dynamics to be 
predominately incoherent over a large reorganisation en- 
ergy rangei^^ii^ 

In the main part of Fig. [5] we plot the transfer rate, 
K+, as a function of reorganisation energy, Ai, calculated 
using Redfield theory (red dashed curve), polaron or 
Forster theory (blue dotted curve) j-A and the variational 
theory with (black circles) and without (black crosses) 
inhomogeneous terms. The inset shows the correspond- 
ing time domain dynamics for Ai = 2 cm"-'^, confirnfing 
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FIG. 6. Energy transfer rate k+ as a function of reorgani- 
sation energy, calculated using Redfield theory (red dashed 
curve), Forster theory (blue dotted curve), and the varia- 
tional theory with (black circles) and without (black crosses) 
bath relaxation terms. The inset shows the corresponding 
time domain dynamics for Ai = 2 cm~^, where the varia- 
tional calculation (black solid curve) is identical both with and 
without inhomogeneous terms. Parameters: V — 2^ cm~*. 



e = 100 cm" 



= 53cm-\ and T = 300 K. 



that the transfer is predominately incoherent, even at 
this small reorganisation energy. 

There are a number of interesting features to be seen 
in the main plot. Firstly, consistent with Fig. |31 the 
transfer rate calculated using Redfield theory fails to de- 
crease significantly at large values of Ai; as expected, 
Redfield theory is inadequate to describe the dynam- 
ics for large system-environment coupling strengthsi^^— 
Secondly, and in contrast, the variational theory (both 
with and without inhomogeneous terms) is capable of in- 
terpolating between the small and large reorganisation 
energy limits, again in accord with the dynamics shown 
in Fig. [21 However, we note that around Ai « 4 cm~* 
the variational theory 'jumps' towards the value pre- 
dicted by Forster (or polaron) theory. This behaviour 
can be attributed to a discontinuous change from non- 
zero to zero renormalised coupling Vr in the variational 
theory, as Ai is increased. For the parameters of Fig. [HI 
above Ai w 4 cm~* the variational minimisation condi- 



tion [Eq. (fT3)) ] suddenly reduces to F(ujk) = 1, and so 
the polaron and variational theories then become equiv- 
alent.*^ Finally, in the present context, perhaps the most 
interesting feature of Fig. [5] is that the addition of in- 
homogeneous terms to the variational theory serves to 
reduce the transfer rate for intermediate values of Ai, 
bringing it closer to that predicted using the hierarchi- 
cal equations of motion or path integral techniques (see 



Refs. |33| and l43l for similar plots using identical parame- 
ters). 

To conclude this section we note that for most param- 
eter regimes the inhomogeneous terms seem to be well 
behaved, and indeed we have just seen how they can in- 
crease the accuracy of the variational method. However, 
we have found that for certain parameters the inhomo- 
geneous terms can be badly behaved, and can even cause 
the variational theory to predict unphysical behaviour 
(e.g. for very large Ai at T = 0), where the behaviour 
is perfectly physical in their absence. We suspect that in 
these regimes the initial bath displacement is too great to 
be described accurately by only the first and second order 
inhomogeneous contributions to the master equation, but 
further work is required to fully assess this conjecture. 



V. SUMMARY 

In summary, we have presented a versatile time-local 
variational master equation technique, which we have 
used to investigate the energy transfer dynamics in a 
donor-acceptor pair. The master equation is constructed 
from a variationally-optimised interaction Hamiltonian, 
obtained by applying a unitary transformation to the full 
Hamiltonian, which results in a perturbative term that 
remains small over a wide range of parameter regimes. 
The formalism also provides a mechanism with which to 
(approximately) include the effects of the dynamic re- 
laxation of the environment. One of the most impor- 
tant aspects of the variational master equation is its flex- 
ibility. We have seen how it naturally reduces to the 
well-known Redfield, polaron, and Forster forms in the 
appropriate limits, and that it can also be used to ex- 
plore energy transfer dynamics under conditions where 
those approaches are known to fail. In particular, we can 
capture coherent dynamics for a donor-acceptor system 
coupled to an Ohmic bath, something not possible in the 
polaron limit. 

The work presented here opens up many potential av- 
enues for future research. For example, although we have 
shown that the qualitative predictions of the variational 
master equation are generally good, it would be interest- 
ing to quantitatively assess its accuracy against known 
numerically-exact benchmarks i^^'^i^'^^ Furthermore, it 
may well be the case that refinements to the variational 
procedure, or even to the form of the unitary transforma- 
tion itself, could increase the accuracy and range of valid- 
ity of the method yet further. An extension to multi-site 
systems would allow a master equation exploration of the 
energy transfer dynamics valid over a range of different 
parameters, for example, in the full FMO complex»2^"— 
Lastly, a study utilising the formalism presented in Ap- 
pendix [B] which extends the theory to include spatially- 
correlated environmental fluctuations, could help to shed 
light on their influence in coherent energy transfer. 
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Appendix A: Master Equation Derivation 

1. Homogeneous terms 

Here we present further details of the derivation of the 
variational master equation used in this work. Upon 
inserting the interaction Hamiltonian into Eqs. (fT6)) 
and (|17p . and returning to the Schrodinger picture, we 
find a master equation of similar form to Eq. (jlSp . 

^ = TTB[IC{t)PxTit)]+TTB[mQxTm, where the 
homogenous contribution can be written 

TrB[IC{t)VxT{t)]^-t[Hs,PTit)] 

ij u> 



with Hs — \t<Jz + Vfiax + i?l. The indices i and j 
{i,j — 1,2,3,4) label system operators Ai = |1)(1|, A2 = 
|2)(2|, A3 = ax, and A4 — ay. In deriving Eq. (|A1|) we 
have utilised the following decomposition of the system 
operators, 

A,,.^ J2 \E){E\A\E'){E'\, (A2) 

E'-E=uj 

where the summation runs over all pairs of eigenval- 
ues of Hs having a fixed energy difference w. Us- 
ing Eq. (jA2p the transformation into the interaction 
picture is easily achieved since we have Aij^{t) — 
exTp[iHst\Aii^exp[—iHst] = A^ ,^ exp[— iwi], while Ai — 
J^ui'^i,'^- Given our form for Hs we find that the u! 
summation runs over the three values uj = 0, ±77, where 
1] — -y/e^ -|- 4yJ as in the main text, and the eigenoper- 



ators are given by Ai ,, 



■ cos 6* sin 61 1 -)(-H I, Ai,o 



■'e\-){-\ +cos^0|+)(+|, A2,, 



cos f sm ( 



|-)(- 



A2,o - cos2 |-)(-| +sin2 6 |+)(+|, A3,,, = cos 20 |-)(+|, 
A3,o = sin20(|+)(+| - |-)(-|), A4,,, - «h)(+|, and 
A4,o — 0. In all cases Ai__,, — A]^, and eigenstates 
of Hs are defined to satisfyi^s |±) = (l/2)(2i? ± t]) |±). 
The angle 9 — (1/2) arctan(2FR/e) characterises the rel- 
ative strength of the renormalised excitonic transfer in- 
teraction to the energy mismatch between the donor and 
acceptor. 
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The time-dependent rates and energy shifts appearing 
in Eq. (|A1I) are given by 'yij{uj,t) = 2Ke[Kij{Lu,t)] and 
Sij{u!, t) = hnlKiji^uj^ i)], respectively, defined in terms of 
the response functions 



Kij{uj,t) 



A,,(T)e*"*dT, 



(A3) 



which depend on the bath correlation functions Aij (r) = 
Tr(B,(T)4(0)pR). Here, we label B, = J2k(9k - 

/fe)(^L + bk,^) for i = 1,2, while B3 = VB^ and 
-B4 — VBy. These correlation functions are given by 
Eqs. dUl) - (HSl) of the main text. 



2. Inhomogeneous terms 

The inhomogeneous terms in the master equation de- 
pend on the quantity Qxt(O). If we wish to consider the 
initial state x(0) = |1)(1| '^ PR, we find 

Sxr(O) = (1 - V){c^ |1)(1| ® prc-G) 

= |1)(1|«>(ptr-Pr)- (A4) 

where pxR = B^ip^B^i is a transformed bath reference 
state. Using Eq. (|17p wc find that the inhomogeneous 
terms take the form 

TTBimQxTiO)] = -z^^r,(i)e^"*[A„ Ai,^] 



{d), 



EE(Hi (^',i)cosc.t-5i;^(c.',t)sin..i) 



X 



[A„A,,^,Ai^-Al,^A],^,] 



EE(Hi (^''Osin..i-|-5l/'(..',i)cosa.i 



x[yl„A,-^,Ai,<. + Al^^A]^^,], (A5) 

in the Schrodinger picture, where Ti{t) — Tr b {Bi{t) p^jx) . 
In a similar way to those appearing in the homo- 
geneous terms, the rates and energy shifts are de- 
fined as -flf{Lo,t) = 2Re[K^^f{iu,t)] and s'^f{uj,t) = 

Ini[_R',- (aj,i)], this time in terms of the response func- 
tions 



Klf(u;,t)^ / A^\T,t)e^-'dT, (A6) 



with A[f{T,t) = TvB{B,{t)B,{t - t)ptr) ~ A,y(r). 

In order to evaluate the quantities Ti{t) and A^ (r, i) 
we utilise the cyclic invariance of the trace and insertions 
of the identity operator, / = _B_|_.ii?_.i, to write 

r. it) = TvB [B, it)pTn) = Trs (B, it)pn) ( A7) 



and 



A|f (r, t) = TTB{B,{t)B,{t - t)pr) - A,,(t), (A8) 



where Bi{t) = B_^iBi{t)B+^i. With reference to Eq. ^, 
we find Biit) - Biit) + c'(i), ^2(0 = Mt), Mt) - 
(F/2)(i?+(i)C+(t) + B^{t)C-{t) - 25), and Bi{t) = 
liV/2){B+{t)C+{t)-B^{t)cS{t)), where 

Ci{t) = 2 f dLu^^^FiLo){l~F(uj))cosLut (A9) 



(AlO) 



and C±(i) — exp[±ii/)(i)], with 



J(c.) 



?/'(t) = 2 / dcj— ^F(w)^sinwi 



We therefore find Ti{t) ^ Ci(t), r2(i) = 0, while r3(t) = 
yR(cos'0(i) — 1), and r4(i) = —ViiSmip{t). In a similar 
way, we find A^f (t, t) = Ci (i)Ci (t - r) while A^j^ (t, t) = 
A21 (r, t) = A22 (t, i) — 0. Those quantities arising from 
^DISPLACED are given by 

V? 



A 



(r, i) = -^ (e^(^) ( cos[^(i) - iPit - t)] ~ l) 

-he-'^(^) ( cos[i^{t) + tPit - t)] - 1) 

-2(cos-0(i)+cosV'(i-r) -2)V 

(All) 



A 



-t)]-1) 



(T,t)=^(e*W(cos[^(i)-^(i-r)] 

-e-'^(^)(cos[7/'(i)-|-^(t-T)]-l)V (A12) 



A 



(r,t) = ^(e^(-)sin[^(f)-^(t-T)] 

e-"^(^) sin[V'(t) + V(i - t)] + 2 sin V^t - r) V (A13) 



and 



A 



(T,t) = ^(-e*Msin[V'(i)-^(t-r)] 
-e-'^(^^ sin[V'(i) -h V'(^ - t)] + 2 sin ^(i) j , 



(A14) 



where (f){t) is defined in Eq. ((22)) . Correlation functions 
of this type coming from the product of -ffoiSPLACED and 
^LINEAR are given by 



Ad) 



A[l>{T,t) = zAc(t) sin^(t - r) + Ti{t)T3{t - r). 



A23^(^,^) = -«Ac(T)sinV'(i-r), 



(A15) 
(A16) 



Ai?(T,t) = -*Ac(T)sinV(i)+ri(t~T)r3(i), (A17) 



A(^2^(T,i)=*Ac(T)sinV'(t), 



(A18) 
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and 

Ad) 



ATJ{T,t) = iAc{T)[cosi;{t - r) - 1] +ri(i)r4(t - r), 

(A19) 



Am (^. t) = -*Ac(t) [cos ?A(i - t) - 1] , 

Ait^(T,t) = -^Ac(T)[cos7^(i) - 1] +ri(t - T)r4(i), 

Ai^^(T,t) = iAc(r)[cosV'W-l]- 

Appendix B: Correlated Environments 



(A20) 



(A21) 
(A22) 



Here we extend the variational formalism presented in 
Section In] to include the effects of spatial correlations in 
the environmental fluctuations. To do so, rather than 
assuming that each site is coupled to a separate bath 
of oscillators, we instead couple each site to a common 
environment. The total Hamiltonian now reads 

ff = ^ ej |A)^.(A| + V{ \XG){GX\ + \GX){XG\ ) 

+ E \^h(^\ T.'^dkA + aljbk) + Hb, (bi) 

j = 1.2 k 

with Hb = X]fc '^fc^I^fci ^^^ ^6 note that we now have 
just one creation operator for each mode, &^. We make 
the spatial separation between the sites explicit with po- 
sition dependent phase factors in the coupling constants. 
We write gkj = 5^6*'° '"^ with rj the position of site j 
(this further assumes each site is coupled to the environ- 
ment with the same strength) . The Hamiltonian describ- 
ing dynamics within the single excitation subspace has a 
similar form as before, 



m 



SUB 



-cr., + V<Tx + Hi 



E m\Y.(9>^A 



ghh), 



(B2) 



i=i,2 



with the Pauli operators defined in the same way as in 
the uncorrected case. 

The variationally-transformed Hamiltonian is once 
again defined as Ht — e'^-ffsuBC"'^, but now we have^ 



G = Eb-)oiE(— ^I-— ^0' 



Jk,_ 



(B3) 



with the variational parameters also containing phase 
factors, fk.j = fk^^^^^ ■ The transformed Hamiltonian 
reads Ht = Hq + Hj, where 

Ho^^ea, + VRa:, + HB + ^Rt, (B4) 

and Hi = i^LiNEAR + -^displaced- Here, 

i^LiNEAR = Yl IJ')(j1 E(5fc - fkKe^'^-^^bl + e-^^-^^bk), 

j k 

(B5) 



and i^DiSPLACED = V{axBx + cTyBy), with B^ = 
{1/2){B+ + B^-2B) and By = {i/2){B+-B_) as before, 
but now 



, ^ l^k ' 



(B6) 



where D{xk) — exp[a;fe6j, — x^bk] is a displacement op- 
erator. Importantly, the renormalisation factor, B = 
Tt{B±pb) now depends on the distance between the two 
sites. 



B = exp r - ^ ^(1 - cos(fc • d)) coth( 



<3"fc ' 



k fe 



(B7) 



where d = ri — r2. The free energy minimisation is 
performed in the same way as in Section [IT] and again 
gives fk = gkF{ujk), but now we find 

F{iOk) = (l + ^-^(l-cos(fc-d))tanh(^)coth(^)) . 
V rjujk / 

(B8) 
The derivation of the variational master equation then 
proceeds in the same way as in the uncorrelated case, and 
results in a form identical to Eqs. (JA1[) and (IA5|) . The dif- 
ference, however, is that the functions Ay(r), A^- (t, i), 
and ri(t), now also have the potential to depend on d. 
We find 

Ai2(r)=A2i(r)= ^3^.(1 -^(a.,))2 

k 

X ( cos(a;fcT — k ■ d) coth(/3a;fc/2) — i sin(ajfcT — k ■ d)) , 

(B9) 

while Aii(t) — A22(t) and are found by setting d = in 
Eq. dm. Once again A33(t) = (y|/2)(e'^(^)-f e-*(^)-2), 
A44(t) = (V^/2)(e*M -e-*(-)), and A34(t) = A43(r) = 
0, but now we have 

0(r) =2 y %F{LUkf{l - cos(fe • d)) 
k^l 
X (cos(a;fcr)coth(/3wfc/2) — isin(ajfeT)). (BIO) 

The final correlation functions appearing in the homoge- 
neous terms are given by Ai4(r) = A42(t) = iAc(r) and 
A24(t) = A4i(t) = -iAc(T) where 

Ac(t) =Vuy] ^F{LOk){l - F{ujk)){l ~ cos(fe • d)) 
k^^ 
X ( sin(a;feT) coth(/3w/c/2) + i cos{uJkT)) , 

(Bll) 

and A3i(t) = Ai3(t) = A32(t) = A23(r) = 0. 

To find the quantities appearing in the inhomogeneous 
terms, we must first find the transformed bath operators, 
B,{t) = B_,,B,{t)B+^i, where B±,i = UkD{fk,i)- We 
now have Bi{t) — Bi{t) + Gi{t), for i = 1,2, where 

„2 

C2{t) = 2y ^F{u:k){l-F{u:k))cos{LOkt-k-d), (B12) 

I ^k 
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and Ci{t) is found by setting d = in Eq. (IB12[) . Lastly, 
Bsit) = {V/2){B+{t)C+(t) + B_(t)C.{t) ~ 2B), and 
Bi{t) = {iV/2){B+{t)C+{t) - B_{t)C-it)), where 



d 



iP{t) = 2 ^ ^FiaJk)"^ ( sin(a;A;t) + sin(wfet - k ■ d)) . 



k '- 



(B13) 

With these transformed operators the functions Ti{t) 
and A^^- (t, i) can be found in a similar manner to Ap- 
pendix ]X1 
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